function [mx, my] = m_p(obj, bx, by, t, x, z0, zc, p, nmax)
%M_P Summary of this function goes here
%   Detailed explanation goes here
    if nargin == 4
        nmax = 10;
    end
    
    if p ~= 0
        [ac, as, fc, fs] = obj.amp_phase(x, z0, zc, p, nmax);
        mx = -by*ac*cos(p*zc*t+fc) + bx*as*sin(p*zc*t+fs);
        my =  bx*ac*cos(p*zc*t+fc) + by*as*sin(p*zc*t+fs);
    else
        [Wp_n, ~, Qp_n, ~] = arrayfun( @(n) obj.WQ_n(x, z0, zc, n, 0), -nmax:nmax);
        W0 = sum(Wp_n);
        Q0 = sum(Qp_n);
        mx=bx*W0-by*Q0;
        my=by*W0+bx*Q0;
    end
end

